#### 0.准备好环境 ####
#Paper：The Relationship between Sleep Disorders and Lower Urinary Tract Symptoms: Results from the NHANES

#下载读取xpt文件的R包“haven”并载入
#install.packages("haven") # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(haven)

#下载用于数据处理包 dplyr, plyr
#install.packages("plyr") # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(plyr)
#install.packages("dplyr") # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(dplyr) # 链接：https://dplyr.tidyverse.org/reference/mutate.html


#下载用于数据快速预览的包 arsenal，使用 tableby 函数
# install.packages("arsenal") # 下载一次后无需再次下载，可以在代码前加入 # 将代码转换为注释
library(arsenal) # 链接：https://cran.r-project.org/web/packages/arsenal/vignettes/tableby.html

#用于加权情况下的分析
library(survey)
setwd("G:/BaiduNetdiskDownload/NHANES/")
#列出下属的 NHANES 数据文件
# list.files()

#### 1. 定位数据模块和变量，获取源数据 ####
##### 1.1 DEMO-人口学数据提取 #####
### 1.1.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/DEMO_D.htm
demo.d <- read_xpt("2005-2006/Demographics/demo_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
demo.e <- read_xpt("2007-2008/Demographics/demo_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(demo.d)
#colnames(demo.e)

# 数据概览使用 tableby, tableby 函数使用的注意事项-1，分类变量可以在前面加上 factor；如果需要分组，则

tab.demo.d <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + RIDSTATR,
                       data=demo.d) 
summary(tab.demo.d, text=TRUE)

tab.demo.e <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + RIDSTATR,
                       data=demo.e) 
summary(tab.demo.e, text=TRUE)

### 1.1.2 提取研究所需要的变量
# 年龄-RIDAGEYR; 性别-RIAGENDR
demo.data.file <- dplyr::bind_rows(list(demo.d, demo.e))
dim(demo.data.file)
demo.data <- demo.data.file[,c("SEQN", "RIDAGEYR", "RIAGENDR")]
# View(demo.data)#查看数据
# 数据概览使用 tableby 

##### 1.2 BMX-身体测量数据提取 #####
### 1.2.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/BMX_D.htm
bmx.d <- read_xpt("2005-2006/Examination/bmx_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
bmx.e <- read_xpt("2007-2008/Examination/bmx_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(bmx.d)
#colnames(bmx.e)

# 数据概览和 NHANES 官方数据做对比
tab.bmx.d <- tableby( ~ BMXBMI,
                      data=bmx.d) 

summary(tab.bmx.d, text=TRUE)

tab.bmx.e <- tableby( ~ BMXBMI,
                      data=bmx.e) 
summary(tab.bmx.e, text=TRUE)

### 1.2.2 提取研究所需要的变量
# BMI-BMXBMI; 
bmx.data.file <- dplyr::bind_rows(list(bmx.d, bmx.e))
dim(bmx.data.file)
bmx.data <- bmx.data.file[,c("SEQN", "BMXBMI")]
# View(bmx.data)#查看数据


##### 1.3 CRP-C反应蛋白数据提取 #####
### 1.3.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/CRP_D.htm
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/CRP_E.htm
crp.d <- read_xpt("2005-2006/Laboratory/crp_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
crp.e <- read_xpt("2007-2008/Laboratory/crp_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

# #colnames(crp.d)
# #colnames(crp.e)

tab.crp.d <- tableby( ~ LBXCRP,
                      data=crp.d) 
summary(tab.crp.d, text=TRUE)

tab.crp.e <- tableby( ~ LBXCRP,
                      data=crp.e) 
summary(tab.crp.e, text=TRUE)

### 1.3.2 提取研究所需要的变量
# C反应蛋白-LBXCRP; 
crp.data.file <- dplyr::bind_rows(list(crp.d, crp.e))
dim(crp.data.file)
crp.data <- crp.data.file[,c("SEQN", "LBXCRP")]
# View(crp.data)#查看数据

##### 1.4 DIQ-糖尿病数据提取 #####
### 1.4.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/DIQ_D.htm
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/DIQ_E.htm
diq.d <- read_xpt("2005-2006/Questionnaire/diq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
diq.e <- read_xpt("2007-2008/Questionnaire/diq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(diq.d)
#colnames(diq.e)

tab.diq.d <- tableby( ~ DIQ010,
                      data=diq.d) 
summary(tab.diq.d, text=TRUE)

tab.diq.e <- tableby( ~ DIQ010,
                      data=diq.e) 
summary(tab.diq.e, text=TRUE)

### 1.4.2 提取研究所需要的变量
# 糖尿病诊断-DIQ010;
diq.data.file <- dplyr::bind_rows(list(diq.d, diq.e))
dim(diq.data.file)
diq.data <- diq.data.file[,c("SEQN", "DIQ010")]
# View(diq.data)#查看数据

##### 1.5 MCQ-合并症数据提取 #####
### 1.5.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/MCQ_D.htm
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/MCQ_E.htm
mcq.d <- read_xpt("2005-2006/Questionnaire/mcq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
mcq.e <- read_xpt("2007-2008/Questionnaire/mcq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(mcq.d)
#colnames(mcq.e)

tab.mcq.d <- tableby( ~ factor(MCQ160B) + factor(MCQ160G) + factor(MCQ160K) + factor(MCQ160C) + factor(MCQ220) + factor(MCQ160F) + MCQ230A + MCQ230B + MCQ230C,
                      data=mcq.d) 
summary(tab.mcq.d, text=TRUE)

tab.mcq.e <- tableby( ~ factor(MCQ160B) + factor(MCQ160G) + factor(MCQ160K) + factor(MCQ160C) + factor(MCQ220) + factor(MCQ160F) + MCQ230A + MCQ230B + MCQ230C,
                      data=mcq.e) 
summary(tab.mcq.e, text=TRUE)

### 1.5.2 提取研究所需要的变量
# 心衰-MCQ160B; 肺气肿-MCQ160G; 慢性支气管炎-MCQ160K; 冠心病-MCQ160C; 癌症或恶性肿瘤-MCQ220; 脑卒中-MCQ160F; 癌症种类（30是前列腺癌）-MCQ230A、MCQ230B、MCQ230C;
mcq.data.file <- dplyr::bind_rows(list(mcq.d, mcq.e))
dim(mcq.data.file)
mcq.data <- mcq.data.file[,c("SEQN", "MCQ160B", "MCQ160G", "MCQ160K", "MCQ160C", "MCQ220", "MCQ160F", "MCQ230A", "MCQ230B", "MCQ230C")]
# View(mcq.data)#查看数据

##### 1.6 BPQ-血压数据提取 #####
### 1.6.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/BPQ_D.htm
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/BPQ_E.htm
bpq.d <- read_xpt("2005-2006/Questionnaire/bpq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
bpq.e <- read_xpt("2007-2008/Questionnaire/bpq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(bpq.d)
#colnames(bpq.e)

tab.bpq.d <- tableby( ~ BPQ020,
                      data=bpq.d) 
summary(tab.bpq.d, text=TRUE)

tab.bpq.e <- tableby( ~ BPQ020,
                      data=bpq.e) 
summary(tab.bpq.e, text=TRUE)

### 1.6.2 提取研究所需要的变量
# 高血压诊断-BPQ020;
bpq.data.file <- dplyr::bind_rows(list(bpq.d, bpq.e))
dim(bpq.data.file)
bpq.data <- bpq.data.file[,c("SEQN", "BPQ020")]
# View(bpq.data)#查看数据

##### 1.7 KIQ-肾脏状况数据提取 #####
### 1.7.1 提取 Component文件
# NHANES官网链接：
#Prostate Conditions：前列腺状况
#https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/KIQ_P_D.htm   
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/KIQ_P_E.htm

# Kidney Conditions - Urology  肾脏状况-泌尿科： 肾脏疾病，尿失禁和夜尿症的个人访谈数据
#https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/KIQ_U_D.htm    
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/KIQ_U_E.htm
#kiq_p 前列腺状况；kiq_u 泌尿科
kiq.p.d <- read_xpt("2005-2006/Questionnaire/kiq_p_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
kiq.p.e <- read_xpt("2007-2008/Questionnaire/kiq_p_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
kiq.u.d <- read_xpt("2005-2006/Questionnaire/kiq_u_d.xpt")
kiq.u.e <- read_xpt("2007-2008/Questionnaire/kiq_u_e.xpt")

#colnames(kiq.p.d)
#colnames(kiq.p.e)
#colnames(kiq.u.d)
#colnames(kiq.u.e)

tab.kiq.p.d <- tableby( ~ factor(KIQ081) + factor(KIQ101),
                        data=kiq.p.d) 
summary(tab.kiq.p.d, text=TRUE)

### 1.7.2 提取研究所需要的变量
#KIQ081: Have trouble starting to urinate（urinary hesitancy）
#KIQ101: After urinating does bladder feel empty（incomplete emptying）
#KIQ005: How often have urinary leakage（urinary frequency）
#KIQ480: How many times urinate in night（nocturia）

urinary.p.data.file <- dplyr::bind_rows(list(kiq.p.d, kiq.p.e))
dim(urinary.p.data.file)
urinary.p.data <- urinary.p.data.file[,c("SEQN", "KIQ081", "KIQ101")]

urinary.u.data.file <- dplyr::bind_rows(list(kiq.u.d, kiq.u.e))
dim(urinary.u.data.file)
urinary.u.data <- urinary.u.data.file[,c("SEQN", "KIQ005", "KIQ480")]
# View(urinary.p.data)#查看数据
# View(urinary.u.data)#查看数据


##### 1.8 SLQ-睡眠状况数据提取 #####
### 1.8.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/SLQ_D.htm
#https://wwwn.cdc.gov/Nchs/Nhanes/2007-2008/SLQ_E.htm
slq.d <- read_xpt("2005-2006/Questionnaire/slq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
slq.e <- read_xpt("2007-2008/Questionnaire/slq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(slq.d)
#colnames(slq.e)

tab.slq.d <- tableby( ~ SLQ060,
                      data=slq.d) 
summary(tab.slq.d, text=TRUE)


### 1.8.2 提取研究所需要的变量
# 睡眠障碍诊断-SLQ060;
slq.data.file <- dplyr::bind_rows(list(slq.d, slq.e))
dim(slq.data.file)
slq.data <- slq.data.file[,c("SEQN", "SLQ060")]
# View(slq.data)#查看数据

##### 1.9 HIQ-保险数据提取 #####
### 1.9.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/HIQ_D.htm
hiq.d <- read_xpt("2005-2006/Questionnaire/hiq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
hiq.e <- read_xpt("2007-2008/Questionnaire/hiq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(hiq.d)
#colnames(hiq.e)


### 1.9.2 提取研究所需要的变量
# 健康保险-HIQ011;  医疗保险-HIQ031A
hiq.data.file <- dplyr::bind_rows(list(hiq.d, hiq.e))
dim(hiq.data.file)
hiq.data <- hiq.data.file[,c("SEQN", "HIQ011", "HIQ031A")]

##### 1.10 OCQ-工作情况 #####
### 1.10.1 提取 Component文件
# NHANES官网链接：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/OCQ_D.htm
ocq.d <- read_xpt("2005-2006/Questionnaire/ocq_d.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"
ocq.e <- read_xpt("2007-2008/Questionnaire/ocq_e.xpt")#参见上述设置默认路径，文件名称后缀不同"d,e"

#colnames(ocq.d)
#colnames(ocq.e)
# 数据概览和 NHANES 官方数据做对比
tab.work.d <- tableby( ~ OCQ265,
                       data=ocq.d) 
summary(tab.work.d, text=TRUE)


### 1.10.2 提取研究所需要的变量
# 健康保险-HIQ011;  医疗保险-HIQ031A
ocq.data.file <- dplyr::bind_rows(list(ocq.d, ocq.e))
dim(ocq.data.file)
ocq.data <- ocq.data.file[,c("SEQN", "OCQ265")]
# View(ocq.data)#查看数据

#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
# 找到权重变量
# DEMO->MEC
#DIQ：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/DIQ_D.htm
#DIQ中没有权重信息：Interview weights should only be used if questionnaire data are analyzed by itself. 
#However if DIQ data are merged with the MEC examination data or laboratory data, the MEC examination weights should be used for analyses. 
#If DIQ data are merged with laboratory sub-sample data, sub-sample weights should be used for analyses.

#MCQ：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/MCQ_D.htm
#MCQ中没有权重信息: A number of additional questions occur in MCQ that are logically related to other questionnaire sections or exam segments. 
#Analysts should refer to these other NHANES sections as part of any data analysis.
#Please refer to the Analytic Guidelines for further details on the use of sample weights and other analytic issues. 
#The Analytic Guidelines are available on the NHANES website.

#BPQ: https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/BPQ_D.htm
#BPQ中没有权重信息

#KIQ_P: https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/KIQ_P_D.htm
#KIQ_P中没有权重信息: These questions are now asked in the Kidney Conditions section of the MEC Interview together with the urinary incontinence and nocturia questions.

#SLQ: https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/SLQ_D.htm
#SLQ中没有权重信息

#CRP：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/CRP_D.htm
#CRP中没有权重信息

#BMX：https://wwwn.cdc.gov/Nchs/Nhanes/2005-2006/BMX_D.htm
#BMX中没有权重信息：The examination sample weights should be used to analyze the body measurement data.

weight.data <- demo.data.file[,c("SEQN", "WTMEC2YR")]
weight.data$WTMEC2YR <- weight.data$WTMEC2YR/2 #权重的计算方式，需要按 Cycle 进行平均

##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.data.file[,c("SEQN", "SDMVPSU", "SDMVSTRA")]

#### 3.合并上述所有数据（把列拼接起来）-Output #### 
output <- plyr::join_all(list(demo.data, weight.data, survey.design.data, 
                              bmx.data, crp.data, diq.data, mcq.data, 
                              bpq.data, urinary.p.data, urinary.u.data, 
                              slq.data, hiq.data, ocq.data), by='SEQN', type='full')

# View(output) # 跑通到这一步，得到这篇文章所用到的全部原始数据
# dim(output) #[1] 20497  

#### 4.变量转换 ####
##### 4.1 重编码BMI为分类变量 ##### 
# 有点类似ifelse 只是省略了很多参数
# max(output$BMXBMI, na.rm = T) 获取最大值排除NA 
bmiBreaks <- c(0, 18.5, 25, 30, max(output$BMXBMI, na.rm = T))
bmiLabels <- c("underweight", "normal weight", "overweight", "obesity")
output$bmiStatus <- cut(output$BMXBMI, breaks = bmiBreaks, labels = bmiLabels, include.lowest = T)


##### 4.2 重编码C反应蛋白为分类变量 #####
output$CRPStatus <- ifelse(output$LBXCRP <= 1, 0, 1)

##### 4.3 构建共病变量 #####
#MCQ160B: 心衰
MCQ160B_transform <- ifelse(output$MCQ160B == 1, 1, 0)
#肺疾病 MCQ160G: 肺气肿；MCQ160K: 慢性支气管炎
MCQ160G_transform <- ifelse(output$MCQ160G == 1 | output$MCQ160K == 1, 1, 0)
#MCQ160C: 冠心病
MCQ160C_transform <- ifelse(output$MCQ160C == 1, 1, 0)
#MCQ220: 癌症或恶性肿瘤
MCQ220_transform <- ifelse(output$MCQ220 == 1, 1, 0)
#BPQ020: 高血压
BPQ020_transform <- ifelse(output$BPQ020 == 1, 1, 0)
# apply apply(X, MARGIN, FUN, ...)
# X: 数组、矩阵、数据框，数据至少是二维的
# MARGIN: 按行计算或按列计算，1 表示按行，2 表示按列
# FUN: 自定义的调用函数
comorbidity_result <- apply(data.frame(MCQ160B_transform, MCQ160G_transform, 
                                       MCQ160C_transform, MCQ220_transform, BPQ020_transform), 1, sum)

output$comorbidity_index <- ifelse(comorbidity_result == 1, 1, ifelse(comorbidity_result < 2, 2, "3 or greater"))
##### 4.4 衍生睡眠状况变量 #####
output$SLQ060[output$SLQ060 %in% c("7", "9")] <- NA
table(output$SLQ060)
output$sleep.disorder <- ifelse(output$SLQ060 == '1', 'Sleep Disorder', 
                                'No Sleep Disorder')


#### 5.根据文章的筛选策略进行筛选 ####
# 原始筛选策略：
#####  5.1年龄>=40，男性(n = 3506)-inAnalysis.index #####

inAnalysis.index <- which(output$RIAGENDR == 1 & output$RIDAGEYR >= 40)
length(inAnalysis.index) # [1] 3506

data.inAnalysis <- output[inAnalysis.index, ]
dim(data.inAnalysis) 

#####  5.2排除患前列腺癌数据(n = 162)-prostate.exclude.index #####
prostate.exclude.index.1 <- which(data.inAnalysis$MCQ230A == 30)
prostate.exclude.index.2 <- which(data.inAnalysis$MCQ230B == 30)
prostate.exclude.index.3 <- which(data.inAnalysis$MCQ230C == 30)
table(output$MCQ230A) #145
table(output$MCQ230B) #16
table(output$MCQ230C) #1

prostate.exclude.index <- c(prostate.exclude.index.1, prostate.exclude.index.2, prostate.exclude.index.3)
length(prostate.exclude.index) # [1] 162

#####  5.3排除患脑卒中数据(n = 209)-data.stroke.exclude #####

stroke.exclude.index <- which(data.inAnalysis$MCQ160F == 1)
length(stroke.exclude.index) # [1] 209

#####  5.4 排除未完成肾脏、尿失禁和夜尿症调查数据(n = 397) #####
#  Paper：questionnaires on kidney/urological, prostate
# 移动检查中心(MEC)访谈的肾脏状况部分数据 & 尿失禁和夜尿症的个人访谈数据

#KIQ081: Have trouble starting to urinate（urinary hesitancy）
#KIQ101: After urinating does bladder feel empty（incomplete emptying）
#KIQ005: How often have urinary leakage（urinary frequency）
#KIQ480: How many times urinate in night（nocturia）

# a <- apply(!is.na(output[,c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')]), 1, any)
# head(a, 20)
# View(output[,c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')])
# a <- is.na(data.inAnalysis[,c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')])
# head(a, 10)
# head(data.inAnalysis[,c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')], 10)
# View(data.inAnalysis[, c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')])
# 任何一个是na的就是拿出来
urinary.exclude.index <- which(apply(is.na(data.inAnalysis[,c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')]), 1, any))
# data.inAnalysis[urinary.exclude.index, c('KIQ081', 'KIQ101', 'KIQ005', 'KIQ480')] # 排除这 4 项有任意一项为 NA 的行

length(urinary.exclude.index) # 397

#####  5.5 排除无睡眠障碍信息的数据 (n = 8) #####
#  Paper：sleep conditions

sleep.exclude.index <- which(is.na(data.inAnalysis$SLQ060))
length(sleep.exclude.index) # 8

##### 5.6 汇总 data.inAnalysis 中要排除的变量(n = 711) #####
exclude.index <- unique(c(prostate.exclude.index, #stroke.exclude.index,
                          urinary.exclude.index, sleep.exclude.index ))
length(exclude.index)  #546
dim(data.inAnalysis) #3506


paper.data <- data.inAnalysis[-exclude.index,]
dim(paper.data) #[1] 2960  


table(paper.data$SLQ060)

